home *** CD-ROM | disk | FTP | other *** search
/ Atari Mega Archive 1 / Atari Mega Archive - Volume 1.iso / language / examples.zoo / numberth / primzahl.lsp < prev    next >
Lisp/Scheme  |  1991-10-22  |  41KB  |  1,038 lines
  1. ; Bruno Haible 10.11.1987-12.12.1987, 25.2.1990
  2. ; Primfaktorzerleger für ganze Zahlen
  3.  
  4. (provide 'primzahl)
  5. ;(in-package 'primzahl)
  6. ;(export '(primteiler pfzv pdivisors divisors pfz prim-von-bis
  7. ;          *Kleinheit* prime probprime notprime factored *pfz-table* 
  8. ;          isprobprime isprime verffz setpfz))
  9. (require 'avl) ; AVL-Bäume
  10.  
  11.  
  12. (setq *print-circle* t) ; wegen der zyklischen Listen in fz-other
  13.  
  14.  
  15. ; (intlist a b) ergibt (a a+1 ... b-1 b), wenn a und b integers sind.
  16. (defun intlist (a b)
  17.   (do ((l nil (cons i l))
  18.        (i b (1- i)))
  19.       ((< i a) l)
  20. ) )
  21.  
  22.  
  23. ; exakter Quotient von Integers, schneller als / :
  24. #-CLISP
  25. (defun exquo (a b)
  26.   (multiple-value-bind (q r) (floor a b)
  27.     (unless (zerop r) (error "Quotient ~S/~S nicht exakt." a b))
  28.     q
  29. ) )
  30.  
  31.  
  32. (require 'jacobi)
  33. ; (jacobi a b) liefert für ganze Zahlen a,b mit ungeradem b>0 das Jacobi-Symbol
  34. ; (a / b). Es ist =0 genau dann, wenn a und b nicht teilerfremd sind.
  35.  
  36.  
  37. ; (pfz n) liefert die Primfaktorzerlegung von n>0, die Primfaktoren von n
  38. ; in aufsteigender Reihenfolge, evtl. mehrfach vertreten.
  39.  
  40. #|
  41. ; vorläufige Version des Primfaktorzerlegers
  42. (defun pfz (n)
  43.   (reverse
  44.     (let ((l nil))         ; l = Liste der bereits gefundenen Faktoren
  45.       (do () ((oddp n))
  46.         (setq n (exquo n 2))
  47.         (push 2 l))
  48.       (let ((p 3))         ; p = Test-teiler (ungerade)
  49.         (loop
  50.           (if (= n 1) (return l))
  51.           (setq s (isqrt n))
  52.           (do () ((or (> p s) (zerop (mod n p))))
  53.                  (incf p 2))
  54.           (if (> p s) (return (cons n l)))
  55.           (push p l)
  56.           (setq n (exquo n p))
  57. ) ) ) ) )
  58. |#
  59.  
  60.  
  61. ; (Miller-Rabin-test N) testet, ob eine gegebene Zahl n>1 wohl eine
  62. ; Primzahl ist.
  63. ; Falls das Ergebnis NIL ist, ist N garantiert zusammengesetzt,
  64. ; Eventuell kommt ein nichttrivialer Teiler von N als zweiter Wert heraus.
  65. ; Falls das Ergebnis T ist, ist N mit hoher Wahrscheinlichkeit prim.
  66. ; (Fehlerwahrscheinlichkeit < 1E-30 )
  67. ; Ist n gerade und >2, so ist das Ergebnis NIL.
  68. (defun miller-rabin-test (n)
  69.   ; Vorausschleife, die bereits 87% aller Zahlen faktorisiert
  70.   (dolist (p '(2 3 5 7 11 13 17 19 23 29 31 37 41 43 47 53 59 61 67))
  71.     (if (zerop (mod n p))
  72.         (return-from miller-rabin-test (if (= n p) t (values nil p)))))
  73.   ; erst jetzt geht's richtig los:
  74.   (let (blist n1 s t1)
  75.     (if (< n 2) (return-from miller-rabin-test nil))
  76.     (setq blist (cond ((< n 2000) '(2))
  77.                       ((< n 1300000) '(2 3))
  78.                       ((< n 25000000) '(2 3 5))
  79.                       ((< n 3200000000) '(2 3 5 7))
  80.                       (t '(2 3 5 7 11 13 17 19 23 29
  81.                            31 37 41 43 47 53 59 61 67 71
  82.                            73 79 83 89 97 101 103 107 109 113
  83.                            127 131 137 139 149 151 157 163 167 173
  84.                            179 181 191 193 197 199 211 223 227 229))))
  85.     ; mit 50 Primzahlen: P(falsch) < (1/4)^50 < 1E-30
  86.     (setq n1 (- n 1))
  87.     (setq s n1  t1 0)
  88.     (do () ((oddp s))
  89.            (setq s (exquo s 2)) (incf t1))
  90.     ; n-1 = s * 2^t1 mit ungeradem s zerlegt.
  91.     (dolist (b blist)
  92.       (if (zerop (mod n b))
  93.         (if (= n b) (return-from miller-rabin-test t)
  94.                     (return-from miller-rabin-test (values nil b))))
  95.       (do ((c (exptmod (mod b n) s n) (sqrmod c n))
  96.            (c1 n1 c)
  97.            (j 0 (1+ j)))
  98.           ((or (= j t1) (= c 1))
  99.            (cond ((not (= c 1)) ; sicher j=t1, aber c/=1
  100.                   (return-from miller-rabin-test nil))
  101.                  ((not (= c1 n1)) ; sicher c=1, c1/=n-1, also j>0
  102.                   ; ggt(c1-1,n) * ggt(c1+1,n) = ggt(c1^2-1,n) = ggt(c-1,n) = n
  103.                   (return-from miller-rabin-test (values nil (gcd n (- c1 1)))))
  104.       )   )) ; Eine Iteration über b beendet
  105.     )   ; blist abgearbeitet
  106.     t   ; n wahrscheinlich prim
  107. ) )
  108.  
  109. ; Hilfsfunktion für Miller-Rabin-Test:
  110. ; Quadrieren modulo n von x mit 0<=x<n
  111. (defun sqrmod (x n)
  112.   (mod (* x x) n))
  113.  
  114. ; Hilfsfunktion für Miller-Rabin-Test:
  115. ; Potenzieren modulo n von a mit 0<=a<n und b>0
  116. (defun exptmod (a b n)
  117.   (let ((c 1))
  118.     (loop    ; a^b*c mod n bleibt invariant
  119.       (if (oddp b) (setq c (mod (* a c) n)))
  120.       (if (= b 1) (return c))
  121.       (setq a (mod (* a a) n) b (floor b 2))
  122. ) ) )
  123.  
  124.  
  125. ; verifizierter Primzahltest für kleine Zahlen n>0
  126. (defun klprimtest (n)
  127.   (cond ((= n 1) nil)
  128.         ((evenp n) (= n 2))
  129.         (t (let ((s (isqrt n)))
  130.              (do ((i 3 (+ i 2)))
  131.                  ((> i s) t)
  132.                (if (zerop (mod n i)) (return nil))
  133. ) )     )  ) )
  134.  
  135.  
  136. ; Primfaktorzerlegung kleiner Zahlen n>0
  137. (defun klpfz (n)
  138.   (let ((l nil))
  139.     (do () ((oddp n)) (setq n (ash n -1)) (push 2 l)) ; Zweier heraus
  140.     (block Rest-prim
  141.       (do ((p 3))
  142.           ((= n 1))
  143.         (loop ; in dieser Schleife wird ein Teiler gesucht
  144.           (cond ((> (* p p) n) (return-from Rest-prim (push n l)))
  145.                 ((zerop (mod n p)) (push p l) (setq n (exquo n p)) (return))
  146.                 (t (incf p 2))
  147.     ) ) ) )
  148.     (nreverse l)
  149. ) )
  150.  
  151.  
  152. ; Datenstruktur: In *pfz-table* steht eine Tabelle von Primfaktorzerlegungs-
  153. ; ergebnissen. Zu jeder Zahl (sie ist der Schlüssel) steht drin:
  154. ; Ob sie Primzahl ist, eine eventuelle Faktorzerlegung bzw. evtl. eine
  155. ; Primitivwurzel.
  156.  
  157. (defparameter *Kleinheit* 10000
  158.   "Nur Zahlen >= *Kleinheit* werden in die PFZ-Tabelle aufgenommen.")
  159. ; *Kleinheit* sollte nur geändert werden, wenn gleichzeitig
  160. ; (setq *pfz-table* nil) ausgeführt wird!
  161. (assert (> *Kleinheit* 1))
  162.  
  163. (deftype fz-class () "Das Ergebnis einer Faktorzerlegung, ganz grob"
  164.   '(member              ; ein Symbol:
  165.             nil        ; (nur kurz nach der Initialisierung)
  166.             prime      ; n ist prim
  167.             probprime  ; n ist wahrscheinlich prim
  168.             notprime   ; n ist garantiert zusammengesetzt, Faktoren unbekannt
  169.             factored   ; n ist zusammengesetzt, schon faktorisiert
  170. )  )
  171.  
  172. (defstruct fz     "Eine Faktorisierung als Ganzes"
  173.   (num 2   :type integer :read-only t)  ; die zu faktorisierende Zahl, >1
  174.   (cl  nil :type fz-class) ; Klassifikation
  175.   (other nil) ; weitere Information:
  176.      ; bei cl = prime : NIL oder eine Primitivwurzel mod n
  177.      ; bei cl = probprime : Eine zyklische Liste von wartenden Funktionen.
  178.      ;                      Jede von ihnen liefert NIL, wenn sie wieder aufge-
  179.      ;                      rufen werden will, und (NIL . factor), falls sie
  180.      ;                      die Zahl als zusammengesetzt nachgewiesen hat, evtl.
  181.      ;                      einen nichttrivialen Faktor gefunden hat, und
  182.      ;                      (T . Prw), falls sie die Zahl als prim nachgewiesen
  183.      ;                      hat, mit evtl. einer Primitivwurzel Prw.
  184.      ;                      (Eventuell zu Beginn die leere Liste.)
  185.      ; bei cl = notprime : Eine zyklische Liste von wartenden Funktionen
  186.      ;                     (Closures), die mitten in der Suche nach einem
  187.      ;                     Faktor von n stecken. Jede von ihnen liefert NIL,
  188.      ;                     wenn sie wieder aufgerufen werden will, und einen
  189.      ;                     nichttrivialen Faktor, wenn sie einen gefunden hat.
  190.      ;                     (Eventuell zu Beginn die leere Liste.)
  191.      ; bei cl = factored : Ein Konstrukt vom Typ fzl, das die Faktorenzerlegung
  192.      ;                     von n widerspiegelt.
  193. )
  194.  
  195. (defstruct fzl "Faktorenzerlegung, Liste"
  196.   (allprimes nil  :type symbol) ; Flag, ob alle Faktoren von factorlist
  197.                                 ; bekanntermaßen Primzahlen sind. 
  198.   (factorlist nil :type list) ; Eine Liste (f1 ... fk) von Faktoren von n,
  199.         ; mit 1 < f1 <= ... <= fk und n = f1 * ... * fk. Die fj sind dabei
  200.         ; (Pointer auf) Konstrukte vom Typ fz, bzw. bei fj<*Kleinheit* ist fj
  201.         ; die Zahl selbst und Primzahl.
  202. )
  203.  
  204.  
  205. ; (fz-key f) ergibt zu einem Element einer Faktorzerlegungsliste die
  206. ; wirkliche Zahl.
  207. (defun fz-key (g)
  208.   (if (integerp g) g (fz-num g))
  209. )
  210.  
  211.  
  212. ; (fz-isprime f) stellt fest, ob ein Element einer Faktorzerlegungsliste
  213. ; prim ist.